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1 Introduction 

Among the different approaches to the 
bioelectromagnetic inverse problem, the current- 
density reconstruction methods (CDR) provide the 
most general solutions. Since the inverse problem 
does not have a unique solution, model assumptions 
have to be taken into account. Multi-channel 
measurements contain not only spatial, but also 
temporal information about the sources, so a 
naturally extension to existing methods leads to 
spatio-temporal model constraints. Spatio-temporal 
CDR’s (stCDR) have been tested in simplified 
volume conductor models, assuming different 
spatial model constraints and a smooth temporal 
activation model. Comparison to existing spatial 
model constraints showed a significant improvement 
of spatial and temporal resolution of the 
reconstructed sources for the spatio-temporal models 
especially in noisy data. 

2 Methods 

2.1 Current density reconstruction (CDR) 

The solution of the bioelectric forward problem is 
given in discretized form by 

m = L*j. (1) 

The dimension of j is 3P, P being the number of 
discretization points in space, N for m, N being the 
number of sensors used in the measurement and 
hence the matrix L, the Leadfield, has the size 
Nx3P. 

The inverse problem is to find the current density for 
a given set of measured values m. Since (1) is 
underdetermined, because there are far more 
discretization points than sensors and the Leadfield 
itself is ill-conditioned, additional constraints on j 
have to be introduced, leading with the most basic 
assumptions about the a-priori state of the current 
density to the minimum norm solution, given by 

j so( (f,) = argmin( 1.• jfi,)-mfi,A |b• j(t,)|p|Y ( 2 ) 

B is the identity-matrix in the case, where minimum 
norm of j without spatial weighting is assumed, and 


is in general of dimension 3Px3P. In the case of 
spatial weighting, B is of diagonal structure and 
contains the weights for each component or location. 
The factor X weights the model assumption versus 
the data. 

The properties of this solution and possible 
variations have been discussed elsewhere [1],[2]. 
Equation (2) fits the source activity to the data and 
to the model specifications, as defined by B, at one 
instance in time. A temporal solution is achieved by 
applying (2) to each measured timeslice of data 
separately. This method, however, has some 
principal disadvantages, e.g., depending on the 
chosen model, the inability to discern close-by 
sources or symmetric source-configurations, which 
can lead to mislocalization of the reconstructed brain 
activity. A further source of reconstruction errors is 
noise in the data, which cannot be suppressed 
completely in the data pre-processing. 

A way to improve the reconstruction results would 
be, to incorporate further model constraints for the 
reconstruction. This bears the danger of 
overspecification of a model, as can be seen in the 
single dipole model, but leads to more robust 
solutions. One can introduce additional constraints 
into the model without loosing generality, by taking 
the time-dependency of the source-activity into 
account. 

2.2 Spatio-temporal current density 
reconstruction (stCDR) 

Each component of j(t) represents several thousand 
neurons being synchronously active at time t. It is, 
from a physiological point of view, reasonable to 
assume, that these ensembles of neurons, that is the 
components of j.(t), do not change abruptly over a 
given interval At of time. Therefore a constraint of 
the temporal behaviour of the current density would 

be — j|(f) —> min From this and the demand, that 

dt ~ 

the reconstructed sources explain the data well at all 
measured timeslices of data, equation (2) can be 
formulated in a time-dependent way: 



J contains all the single timeslice solutions j, and M 
likewise the measured data for all timeslices, 

j r =(j(^i),j(^X---,j(^)) rand 

S is the number of samples or timeslices. The 
dimensions of the matrices L T , B T and K, the 
forward model operator, the spatial model operator 
and the temporal model operator are the same as in 
equation (2) multiplied by the number of timeslices 
T. The structure of equation (3) is the same as that 
of (2), the problem has merely been blown up by a 
factor T, so all numerical procedures that solve (2) 
can also, in principle, be applied to equation (3). 
Since the head-model, from which L is generated, is 
time-independent, L T is of block-diagonal structure 
and each block contains L. The same applies to B T 
,if no spatio-temporal correlation or coupling 
between the components of J is assumed. Here the 
blocks of B t contain the desired spatial model 
operator B. 

K is given by 
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if minimal temporal change in j is desired and is the 

discrete formulation of — j (t ), as proposed by M. 
dt~ 

Wagner [2], 

2.3 Solution of the linear problem 

If p=2, the solution of (3) reduces to the solution of 
linear system. The corresponding normal equation of 
(3) is 

(L; *L +Aj;*B +lp .K).j r =L;( 4 ) 

The solution of this can be given in matrix 
formulation by 

S»J + J«T = L «M , 

S = L' •L + A S B' *B,T = A r W' *W 


j = 0(^1 )p2 )■■■ jOs)), M = (m(^) m (t 2 )■■■ m(t s )) 

Equation (5) is of the form the continuous Sylvester- 
equation and can be solved by direct means in the 
case of a small number of influence points P. In 
realistic head-models, the number of influence 
points is large. For a 5mm discretization of the brain 
by a regular grid, there are P~12000 possible source 
locations. Eq. (5) is then of order 36000 and not 
solvable within reasonable amounts of memory and 
compuation time. In this case, an iterative solution 
of the problem can be applied, as proposed in [3], 

The normal equation of the problem has a block 
structure 



with blocks 

A fi =L‘ •L + A S B‘ »B +A t (W* • W) a »l 3p , 

A| = A t (W‘ • W)ij •ljp 

The iterationscheme is then given by 
j^^A^L'.mi-gA .f) (5) 

Since the A., have a sparse structure, (5) can be 
solved with the conjugate-gradient-method, even in 
the case of large P. 

An effcient formulation of the problem can be 
found, if, like the spatial coupling in LORETA [4], 
the temporal coupling is applied to the weighted 
sources. 

The normal equation for the problem is then 
(L;.L +Ajf «B + 4 jf r *Kf •K.B ).j =L; .m / ( 6 ) 
and can be rewritten as 

(£ =£ ( ? )’ 

L =L B 1 , j = B j 

T =T=T ’ ±T T J T 

The matrices L t and K commutate and thus the 
solution of (7) is equal to 


( 5 ) 






(L«l‘ +A s l 7v )»U + A r U*(W' • W) = M, 

J = L* «U 

Equation (8) is only of order N and can be solved 
very fast. However, the solution found by (8) is 
different from the original problem. 

3. Simulation 

For testing the properties of the temporal current 
density reconstruction a simple volume conductor 
model was chosen for various reasons: 

• The computational demands with respect to 
memory usage and computation time for stCDR 
are extremely high. In a simple model with a 
small number of discretization points for j, a 
large number of calculations can be carried out 
in a reasonable amount of time. 

• The properties of the simple volume conductor 
are well understood. 

• In order to prove the principle advantages of 
stCDR, the construction of an appropriate depth¬ 
weighting, which is necessary for a source space 
model in a volume, can be neglected in a planar 
source space model. The properties of the 
simplified forward model with respect to 
inversion , especially the bad condition of the 
matrix of the foward-model operator L, are 
similar to those in realistic models. 

3.1 Infinite homogenous volume conductor 

The simple model set up consisted of a 10 by 10 
grid, with a length of 10 arbitrary units per side 
centered at (5.5, 5.5,0.0). Nine EEG-sensors were 
placed in a square planar array with centre at 
x=5.5,y=5.5 and z=2 above the grid (Fig.: 1). The 
forward model was calculated by 


where r ( =is the position of the i-th sensor, r, the 
position of the j-th gridpoint, i.e., we number the 
sensor positions by i and the gridpoints by j. 

This formula applies to an infinite homogenous 
volume conductor. 

Two point sources were placed on the grid and a 
forward solution using (9) was calculated. To each 
point source a gaussian activation was assigned with 
a width of w=2 and center at tp ea k: 


q(0 = q n exp(- - —) 

- - u w 

The simulation was calculated for 20 timeslices. 



3.2 Forward simulations 

The forward model was computed for the following 
conditions: 

• Noise with 4 variations (0, 10, 20, 30% of the 
maximum signal magnitude). The average 
signal-to-noise ratio per channel resulting from 
these noise-levels ranges from 2.0 (10% noise) 
to 1.4 (30% noise). However, since the noise- 
amplitude is derived from the maximum signal 
amplitude, the SNR depends strongly on the 
individual source configuration. Therefore, the 
percentage of the maximum amplitude was 
taken as a measure for the noise instead of the 
resulting SNR. 

• Location of the two sources on the grid with 5 
variations (9, 7, 5, 3, 2 points on the grid 
between sources), 

• Time shift between the two sources with 3 
variations, 

• Orientation of the sources with 3 variations. 

The simulation procedure is shown in figure 2. First 
the data is generated by the original source 
configuration, then noise is added and the 
reconstrucion is done for all timeslices. The 
reconstructed source timeseries are then generated 
form the spatial reconstruction. 







4 Results 


Goodness of temporal reconstruction 


All 180 configurations were tested for: 

• Spatial localization error 

• Goodness of temporal reconstruction 



standard solution coupled solution 



Reconstructed source-time-series 
Figure 2: Steps o/the simulation.On the leftside, the 
results for the 12 solution are shown, on the right 
side the results for 12-coupled solution. The 
contourplots show the reconstructed current density 
for each timeslice. The original source-positions are 
marked by white squares at the times of peak 
activity. The reconstructed timeseries for both 
sources are shown below the reconstruction. 


4.1 Spatial localization error 

In order to determine the spatial localization error of 
a trial, the local maxima of the strength of the 
current density were determined for each timeslice 
separately, beginning at the peak of the first original 
source and ending at the peak of the last original 
source. Each local maximum was interpreted as a 
reconstructed source. The spatial distance of the 
reconstructed sources to the original sources was 
calculated and the two reconstructed sources with 
the shortest distance to the original sources were 
assumed to be the correct reconstructed sources for 
that timeslice. 

The average distance for all sources over all 
timeslices was calculated and termed the spatial 
localization error for that trial: 


A second measure of the goodness of the 
reconstruction is the correlation between the 
reconstructed source-time-series and the original 
source-time-series. The reconstructed source-time- 
series were derived from the reconstructed sources 
at the respective peak times. For the first original 
source, 

Figure 3: Reconstruction errors and goodness of 
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temporal reconstruction for the coupled (dashed) 
and the uncoupled (solid) solution for different 
noiselevels 

the position of the reconstructed source was 
determined at t pea ki and the time-series of the strength 
at that position was taken and the correlation with 
the time-series of the original source was calculated. 

The same procedure was applied to the second 
source. The mean of the two correlation coefficients 
was termed the goodness of temporal reconstruction 
for that trial. 

Figure 4: Reconstruction errors and goodness of 
average localization error goodness of temporal reconstruction 
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temporal reconsunstruction for the coupled (dashed) 
and the uncoupled (solid) solution for different 
source-distances. 

4.3 Conclusion 

Figures 3-6 show the localisation error and goodness 
of temporal reconstruction for the tested 
configurations. 
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Figure 5: Reconstruction errors and goodness of 
temporal reconsunstruction for the coupled (dashed) 
and the uncoupled (solid) solution for differents 
timeshifts. 


additional information, which is provided by the 
time-dependent model constraint (smoothness in 
time). This general advantage of spatio-temporal 
current density reconstruction (stCDR) was 
shown to be most prominent in the case of noisy 
data due to the low-pass filter properties of the 
algorithm, but also in the ability to separate 
sources 

This leads to the conclusion, that, with 
appropriate depth weighting, the stCDR could 
also improve reconstruction in realistic volume 
conductor models. 
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Figure 6: Reconstruction errors and goodness of 
temporal reconsunstruction for the coupled (dashed) 
and the uncoupled (solid) solution for differents 
diretions 


It could be shown in simulations, that the introduction 
of temporal constraints to existing CDRs leads to 
significant improvements in spatial and temporal 
accuracy. This improvement is gained by the 
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